setwd("E:\\Rbigwork\\")
ds <- read.table(file = "GSE42546_CleanedRawCounts.txt",header = T,row.names = 1)

###数据特征可视化###
#柱状图查看总体数据分布#
ds2 <- as.matrix(log2(ds+1))  #标准化
p <- hist(ds2,main = "Histogram of GSE42546",xlab = "log2(data+1)",col = cm.colors(16))
png("histogram.png")
p
dev.off()

#箱线图查看四个分组的数据特征#
library('reshape2')
ds1 <- melt(ds)
colnames(ds1)[1] <- "sample"
ds1[,2] <- log2(ds1[,2]+1)  #标准化
group <- gsub("\\d","",ds1[,1])
table(group)
group1 <- factor(group,labels = c("bipolar disorder",  #B表示bipolar disorder
                                  "control",           #C表示control
                                  "major depression",  #D表示major depression
                                  "schizophrenia"))    #S表示schizophrenia
table(group1) 
ds1$group <- group1  #分组信息
library(ggplot2)
library(ggsci)
p1 <- ggplot(ds1, aes(x=group,y=value,fill=group)) +   #根据group分组绘制箱型图
  geom_boxplot()+
  stat_boxplot(geom = "errorbar",
               lwd=0.5,
               width=0.2)+
  theme_minimal()+
  ggtitle("boxplot of GSE42546")+
  scale_fill_nejm()+
  theme(legend.position ="none")
pdf("boxplot.pdf")
p1
dev.off()

#样本聚类信息#
d1<-dist(t(ds)) #欧氏距离
fit1<-hclust(d1,method = "average")
p2 <- plot(fit1,hang = -1,cex=.8,main ="Sample Hclust") #树状图
print(p2)
p6 <- heatmap(as.matrix(d1))  #热图
print(p6)

d2<-dist(t(ds),method = "manhattan") #曼哈顿距离
fit2<-hclust(d2,method = "average") 
p5 <- plot(fit2,hang = -1,cex=.8,main ="Sample Hclust") #树状图
print(p5)
p7 <- heatmap(as.matrix(d2))  #热图
print(p7)

#主成分分析#
library('ggfortify')
group <- substr(colnames(ds),1,1)
table(group)
group <- factor(group,labels = c("bipolar disorder",  
                                  "control",           
                                  "major depression",  
                                  "schizophrenia"))    
table(group)
ds3 <- as.data.frame(t(ds2))
ds3$group <- group
p3 <- autoplot(prcomp(ds3[,1:(ncol(ds3)-1)] ), data=ds3,colour = 'group',
               frame = T, frame.type='norm',) + 
  theme_classic()
print(p3)


###差异表达分析###
library(DESeq2)
colData <- data.frame(row.names=colnames(ds), group)  #建立样本信息
dds <- DESeqDataSetFromMatrix(countData = round(ds),colData = colData,design = ~group)
dds$group<- relevel(dds$group, ref = "control")   #指定control作为对照组
dds2 <- DESeq(dds)
res <- results(dds2)
summary(res)   
res1 <- data.frame(res) #将res转化为表格格式
write.csv(res1,'DESeq2_result.csv')

res1 <- res1[,c("log2FoldChange","pvalue")]
# 获取pvalue<0.05，差异倍数>1.5或者<-1.5的显著差异表达基因
res1_up<- res1[which(res1$log2FoldChange >= log2(1.5) & res1$pvalue < 0.05),]      # 表达量显著上升的基因,共47个
res1_down<- res1[which(res1$log2FoldChange <= -log2(1.5) & res1$pvalue < 0.05),]    # 表达量显著下降的基因，共210个

#可视化：火山图#
library(ggplot2)
library(ggrepel)
fccutoff=1.5
pcutoff=0.05

#显著上调、下调基因分组
res1$log10pvalue <- -1*log10(res1[,"pvalue"])
res1$group <- "nosig"
pos1 <- res1$log2FoldChange >= log2(fccutoff) & res1$pvalue < fccutoff
res1$group[pos1] <- "up"
pos2 <- res1$log2FoldChange <= -log2(fccutoff) & res1$pvalue < fccutoff
res1$group[pos2] <- "down"

#分别筛选差异显著前6个的上下调基因，并合并两组数值
library(dplyr) # 用于数据处理
library(gt) # 制作表格
top_20 <- bind_rows(
  res1 %>%
    filter(group == 'up') %>%
    arrange(pvalue, desc(abs(log2FoldChange))) %>%
    head(6),
  res1 %>%
    filter(group == 'down') %>%
    arrange(pvalue, desc(abs(log2FoldChange))) %>%
    head(6)
)
top_20 %>% gt()  #将数据制成表

#作图
p4 <- ggplot(res1,aes(x=log2FoldChange,y=log10pvalue)) +
  geom_point(aes(color = group),alpha=0.65, size=2)+
  scale_color_manual(values = c("#546de5", "#d2dae2","#ff4757"))+
  geom_hline(yintercept = -log10(pcutoff),lty=2) + 
  geom_vline(xintercept = c(log2(fccutoff),-log2(fccutoff)),lty=2)+
  labs(y="-log10Pvalue")+
  ggtitle("Volcano plot")+
  theme_bw()
print(p4)

#为差异显著前6个上下调基因添加标记
p4 <- p4 +
  geom_label_repel(data = top_20,
                   aes(log2FoldChange, log10pvalue, label = rownames(top_20)),
                   size = 3, fill="#CCFFFF")

print(p4)

###富集分析###
library(org.Hs.eg.db)
library(clusterProfiler)
#将基因名称进行格式转化#
#显著下调基因
ids1<- bitr(rownames(res1_down), 
           fromType = "SYMBOL",   #输入的格式为SYMBOL
           toType = c("GENENAME","ENTREZID"),   #转化为ENTREZID格式
           OrgDb = 'org.Hs.eg.db')  
#显著上调基因
ids2<- bitr(rownames(res1_up), 
            fromType = "SYMBOL", 
            toType = c("GENENAME","ENTREZID"), 
            OrgDb = 'org.Hs.eg.db')
#背景基因
univ <- bitr(rownames(ds),
             fromType = "SYMBOL", 
             toType = c("GENENAME","ENTREZID"), 
             OrgDb = 'org.Hs.eg.db')

#GO富集分析#
egodown <- enrichGO(ids1$ENTREZID, "org.Hs.eg.db", 
                  pAdjustMethod = "none",  #富集结果不明显，进行修改
                  qvalueCutoff = 1,
                  keyType = "ENTREZID", #设置ID类型为ENTREZID
                  ont = "ALL" ,#选择注释类型为全部
                  universe = univ$ENTREZID  #设置背景基因
)
dotplot(egodown,showCategory=20)  


egoup <- enrichGO(ids2$ENTREZID, "org.Hs.eg.db", 
                    pAdjustMethod = "none",  #富集结果不明显，进行修改
                    qvalueCutoff = 1,
                    keyType = "ENTREZID", #设置ID类型为ENTREZID
                    ont = "ALL" ,#选择注释类型为全部
                    universe = univ$ENTREZID  #设置背景基因
)                    
dotplot(egoup,showCategory=20)

#KEGG富集分析#
keggdown <- enrichKEGG(gene = ids1$ENTREZID, 
                   organism = "hsa",
                   keyType = "kegg", 
                   pAdjustMethod = "none",  #富集结果不明显，进行修改
                   qvalueCutoff = 1,
                   universe = univ$ENTREZID)
barplot(keggdown,showCategory=20,drop=T)

keggup <- enrichKEGG(gene = ids2$ENTREZID, 
                       organism = "hsa",
                       keyType = "kegg", 
                       pAdjustMethod = "none",  #富集结果不明显，进行修改
                       qvalueCutoff = 1,
                       universe = univ$ENTREZID)
barplot(keggup,showCategory=20,drop=T)

#拓展：GSEA富集分析#
#提取降序排列后的log2FoldChange,且名称为ENTREZID
res2 <- data.frame(SYMBOL=rownames(res1),log2FC=res1[,"log2FoldChange"])
genes <- merge(res2,univ,by='SYMBOL')
genelist <- genes$log2FC
names(genelist) = genes$ENTREZID
genelist <- sort(genelist, decreasing = TRUE)
head(genelist)

#KEGG的GSEA分析
gseaKEGG <- gseKEGG(geneList     = genelist,
                    organism     = 'hsa',
                    nPerm        = 1000,
                    pvalueCutoff = 0.05,
                    pAdjustMethod = "none",
                    verbose      = FALSE)
#作图分别查看激活、抑制的通路
dotplot(gseaKEGG,showCategory=12,split=".sign")+facet_grid(~.sign)
#把富集的结果转换成data.frame
gseaKEGG_results <- gseaKEGG@result
#被激活通路：MAPK signaling pathway
library(enrichplot)
pathway.id = "hsa04010" #MAPK signaling pathway的ID
gseaplot2(gseaKEGG, 
          color = "red",
          geneSetID = pathway.id,
          pvalue_table = T)
#被抑制通路：Endometrial cancer
pathway.id = "hsa05213" #Endometrial cancer的ID
gseaplot2(gseaKEGG, 
          color = "red",
          geneSetID = pathway.id,
          pvalue_table = T)
#所有通路，不显示p值
gseaplot2(gseaKEGG,1:9,color="red") 

#pathview查看MAPK signaling pathway基因变化情况
library(pathview)
pathway.id = "hsa04010"
pv.out <- pathview(gene.data  = genelist,
                   pathway.id = pathway.id,
                   species    = "hsa")

